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(54) Method and apparatus for providing real-time calculation and display of tissue deformation 
in ultrasound imaging 



(57) An ultrasound system and method for calcula- 
tion and display of tissue deformation parameters are 
disclosed. An ultrasound acquisition technique that 
allows a high frame rate in tissue velocity imaging or 
strain rate imaging is employed. With this acquisition 
technique the same ultrasound pulses are used for the 
tissue image and the Doppler based image. A sliding 
window technique (160, 161, 162, 163 and 164) is used 
for processing. The tissue deformation parameter strain 
is also determined by an accumulation of strain rate 
estimates for consecutive frames over an interval. The 
interval may be a triggered interval generated by, for 
example, an R-wave in an ECG trace. The strain calcu- 
lation (150) may be improved by moving the sample vol- 
ume from which the strain rate is accumulated from 
frame-to-frame according to the relative displacement of 
the tissue within the original sample volume. The rela- 
tive displacement of the tissue is determined by the 
instantaneous tissue velocity of the sample volume. An 
estimation of strain rate (150) based upon a spatial 
derivative of tissue velocity is improved by adaptively 
varying the spatial offset, dr. The spatial offset, dr, can 
be maximized to cover the entire tissue segment (e.g., 
heart wall width) while still keeping both of the sample 
volumes at each end of the offset within the tissue seg- 
ment. This may be accomplished by determining 
whether various parameters (e.g., grayscale value, 



absolute power estimate, magnitude of the auto-corre- 
lation function with unity temporal lag and/or magnitude 
of strain correlation) of the sample volumes within in the 
spatial offset are above a given threshold. Strain rate 
may be estimated (150) using a generalized strain rate 
estimator that is based on a weighted sum of two-sam- 
ple strain rate estimators with different spatial offsets. 
The weights are proportional to the magnitude of the 
strain rate correlation estimate for each spatial offset, 
and thus reduce the effect of noisy, i.e. poorly corre- 
lated, samples. An improved signal correlation estima- 
tor that uses a spatial lag in addition to the usual 
temporal lag is disclosed. The spatial lag is found from 
the tissue velocity. The improved signal correlation esti- 
mator can be utilized both in the estimation of strain rate 
and tissue velocity. Tissue velocity may be estimated in 
a manner that reduces aliasing while maintaining spatial 
resolution. Three copies of a received ultrasound signal 
are bandpass filtered at three center frequencies. The 
middle of the three center frequencies is centered at the 
second harmonic of the ultrasound signal. A reference 
tissue velocity is estimated from the two signals filtered 
at the outside center frequencies. The reference tissue 
velocity is used to choose a tissue velocity from a 
number of tissue velocities estimated from the signal 
centered at the second harmonic. A method to estimate 
(150) the strain rate in any direction, not necessarily 
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along the ultrasound beam (144), based on tissue velocity data from a small region of interest around a sample volume 
is disclosed. Quantitative tissue deformation parameters, such as tissue velocity, tissue velocity integrals, strain rate 
and/or strain, may be presented (152) as functions of time and/or spatial position for applications such as stress echo. 
For example, strain rate or strain values for three different stress levels may be plotted together with respect to time over 
a cardiac cycle. Parameters which are derived from strain rate or strain velocity, such as peak systolic wall thickening 
percentage, may be plotted with respect to various stress levels. 
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Description 



[0001] The present invention relates to diagnostic ultrasound systems which measure and image anatomical struc- 
tures, and their movement. More particularly, the present invention relates to a signal processing method and apparatus 
5 for calculation and display of tissue deformation to be used in ultrasonic imaging systems. 

[0002] Recently, within the field of ultrasound imaging, physicians have become interested in using tissue deforma- 
tion properties, such as tissue strain and strain velocity for clinical measurements. 

[0003] The term "strain" refers to a characteristic of material being examined. For example, the strain associated 
with muscle tissue corresponds to a ratio of the change in muscle tissue length during a prescribed time interval to the 

w muscle tissue's initial length. In ultrasound imaging, the rate of change of strain (e.g., strain rate, strain velocity, etc.) 
may be visually presented to a physician as a colorized 2-dimensional image, where variations in color correspond to 
different strain velocities. It has become apparent that the viability of a segment of the muscle is related to the amount 
of muscle strain and temporal behavior of the strain that is performed by, or is imposed on the muscle segment. Also, it 
has been determined that malignant tumors may be detected based on their resistance to compression. 

15 [0004] One application of real-time strain velocity imaging is in cardiology. The strain velocity gives a direct and 
quantitative measure for the ability of the myocardium to contract and relax. By imaging along the myocard from an api- 
cal view, the local strain velocity component along the long axis of the heart can be measured. Measuring the local 
strain velocity component gives information about the local shortening and lengthening of the heart wall. By imaging 
from the parasternal view, the strain velocity component perpendicular to the heart wall can be found. Finding the strain 

20 velocity component perpendicular to the heart wall gives information about the local thickening of the muscle. Wall thick- 
ening measured with M-mode or from the 2D image is a commonly used measure for muscle viability. With strain veloc- 
ity imaging, a direct measure for this thickening is available. The strain velocity images can potentially add to the 
diagnosis of a number of cardiac disorders. 

[0005] Another application of strain velocity imaging is in heart transplants. Velocity variations inside the myocar- 
25 dium are important for the diagnosis of rejection after heart transplantation. The strain velocity images give a direct dis- 
play of these velocity variations. 

[0006] Another application of strain velocity imaging is in noninvasive electrophysiology. The preferred embodiment 
describes techniques to image the local contraction/relaxation contributions with a high spatial and temporal resolution. 
Local contraction/relaxation information can be used to accurately determine the localization of, for example, where the 
30 mechanical movement in the heart chambers is activated based on a cross section just below the AV-plane. Further- 
more, aberrant conduction pathways (Wolf-Parkinson-White) from the atrium to the ventricle can be localized for later 
ablation. Even the depth inside myocard of these paths can be better localized with this invention in order to determine 
if the patient should be treated with catheter techniques or surgical techniques. 

[0007] Another application of strain velocity imaging is in measuring cardiac wall thickening. A well established 
35 methodology in cardiac diagnosis is to acquire a M-Mode image and to measure the wall thickening of myocardium dur- 
ing systole. The preferred embodiment provides techniques to take this wall thickening information and measure it in 
real-time with a high precision in both the spatial and temporal domain. The high diagnostic relevance of the current wall 
thickening measurements indicates that the imaging modality described in this invention contains highly relevant infor- 
mation for cardiac diagnosis 

40 [0008] To understand strain velocity or strain rate in more detail, it is assumed that an object of initial length L 0 may 
be stretched or compressed or itself lengthens or contracts to a different length L The one-dimensional strain, defined 
as 



represents a dimensionless description of the change. If the length L is considered to be a function of time, the temporal 
derivative of the strain, the strain velocity, can be found using the equation: 



55 [0009] If the velocity, v of every point in the object is known, an equivalent definition of the strain velocity is: 
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[0010] These equations also provide a useful description of the deformation of the object. In Eq. 3, r is the spatial 
direction of the stretching or compression. The relation between Eq. 2 and Eq. 3 can be seen if the length L is defined 
as L(f) = r 2 (t)-r 1 (t) , and L 0 =L(f 0 ) , where ^ is the distance to one end of the object, and r 2 is the distance to the 
other, and letting t -> to, and letting r1 -> r2. As illustrated in Eq. 3, the strain velocity is in fact the spatial gradient of the 
5 velocity. The strain velocity thus measures the rate of the deformation of the object. If the strain velocity is zero, the 
shape of the object is not changing. If the strain velocity is positive, the length of the object is increasing, and if the strain 
velocity is negative, the length of the object is decreasing. Strain velocity is also known as rate-of-deformation, stretch- 
ing, strain rate or velocity strain. 

[0011] Strain imaging is currently an established research area in ultrasonic imaging. The degree of deformation in 

10 imaged structure may be estimated by correlation of 2D images obtained before and after a pressure increase. One dis- 
advantage of estimating image deformation based on correlation of images is that the instantaneous value of the strain 
is not calculated nor displayed in real-time. The lack of a real-time capability is an important clinical disadvantage. For 
example, if strain imaging could be performed in real-time, strain imaging could be applied more effectively in cardiac 
ultrasound or could be used as an interactive inspection modality where anomalies in tissue compressibility can be vis- 

15 ualized in real-time according to the pressure gradients that are applied to the imaged structures. 

[0012] A method of position tracking has been proposed to estimate the local myocardial strain velocity based on 
radio frequency (RF) M-Mode acquisitions. The position tracking method is described in H. Kanai, H, Hasegawa, N. 
Chubachi, Y. Koiwa, and M. Tanaka, "Noninvasive evaluation of local myocardial thickening and its color-coded imag- 
ing," IEEE Trans, on Ultrasonics, Ferroelectrics and Frequency Control, vol. 44, pp. 752-768, 1997. However, the 

20 method described in the Kanai et al. article has the disadvantages of poor temporal resolution and high computational 
cost, which render real-time imaging difficult and costly. Furthermore, the method described in the Kanai et al. article is 
a manual M-mode technique, not well suited to form the basis for real-time two-dimensional strain images. Also, the 
strain velocity is a derivative of a velocity estimate and is therefore very noise sensitive. The fundamental velocity alias- 
ing problem that is inherent in tissue velocity imaging makes noise difficult to overcome because aliasing prevents the 

25 pulse repetition frequency from being set at a low enough rate to allow a large observation time. If the observation time 
could be increased, the noise robustness of the strain velocity images could be significantly improved. 
[0013] Certain of the above identified difficulties are addressed and overcome according to the teachings of United 
States Patent Application No. 09/167,896, filed October 7, 1998 and entitled "A METHOD AND APPARATUS FOR 
PROVIDING REAL-TIME CALCULATION AND DISPLAY OF STRAIN IN ULTRASOUND IMAGING," However, is an 

30 object of the present invention to supplement and/or improve upon such teachings. Certain additional difficulties and 
shortcomings of the prior art are described below. 

[0014] To achieve high frame rate in color Doppler applications, two previously known techniques are commonly 
used: multi line acquisition (MLA) and interleaving. These techniques make it possible to acquire more data than in a 
basic mode, where the scanner after one pulse is received waits the specified pulse repetition time (T) before firing the 
35 next pulse in the same direction. The time to acquire a frame of Doppler data in the basic mode is: 



where N is the number of pulses in each direction and N b is the number of beams in the image. A relatively small extra 
delay related to the change in setup of the transmitter and beamformer is ignored to simplify the discussion. 
45 [0015] In the MLA method, a broad beam is transmitted. When receiving the echo, the signals from all the trans- 
ducer elements are processed in parallel in two or more beamformers. Each beamformer time delays the element sig- 
nals differently to generate different receive beams. This way, two or more beams can be acquired during the time for 
one pulse-echo cycle, and the frame rate can be increased correspondingly. Using MLA, the time to acquire a frame of 
Doppler data is 



t D0 =N b NT, 



(4) 
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where N MLA is the number of beams that are processed in parallel. 

[0016] In the interleaving technique, the waiting time T from one pulse to the next in the same direction is utilized 
to send pulses in other directions, as illustrated in Fig. 1. There is however a minimum waiting time T 0 where no other 
pulses can be fired in any direction. This is given by the time for the pulse to travel to the maximum depth and back: T 0 
5 > 2dlc. The number of directions that pulses are fired during the time T is called the interleave group size, N int . This 
obviously has to be an integer number, and T = N int T 0 . Using interleaving, the time to acquire a frame of Doppler data 
becomes: 



10 
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[0017] Fig. 1 illustrates the pulse order and beam directions in the interleaving method for three different interleave 
group sizes N int . The number of beams, N b , equals 8 and packet size, N, equals 2 in the example of Fig. 1. For inter- 
leave pattern 100, the interleave group size, N int , equals 8, for interleave pattern 110, N jnt equals 4 and for interleave 
pattern 120, N int equals 1. 

[0018] A typical scanning procedure for a tissue Doppler application is illustrated in Fig. 2. In the example of Fig. 2, 
the packet size N equals 3 and the interleave group size Nint equals Nb. T is the pulse repetition time, t T and t D are the 
times needed to acquire a tissue frame and a Doppler frame respectively, and t F is the total acquisition time for one tis- 
sue Doppler frame. A tissue frame 130 is first captured, using high beam density. The PRF used for tissue Doppler is 
usually so low that only one interleave group is necessary. So N Doppler subframes 132, 134 and 136 are captured 
separately, usually using fewer beams than in the tissue frame. The velocity is calculated from the N subframes 132, 
134 and 136, color coded and then mapped onto the tissue frame. The time to acquire a tissue Doppler frame then 
becomes: 



N 



NT, 



MA 



((7) 



35 



where t T is the time required to acquire the tissue frame. It is, thus, apparent that the maximum frame rate is limited by 
the above described ultrasound data acquisition schemes. 

[0019] It is known that tissue velocity can be estimated using either the first or the second harmonic component of 
40 the ultrasound signal. Using the second harmonic component (octave imaging) has been reported to produce an 
improvement in image quality in gray scale images, and the same improvement can be expected in tissue Doppler. 
There is however a disadvantage in that the Nyquist limit is halved when using the second harmonic instead of the fun- 
damental component. Using a low PRF is also preferable, since the phase amplitude of the complex signal is increased 
compared to the noise, resulting in a lower variance in the velocity estimate. A disadvantage of using a low PRF is that 
45 the Nyquist limit is further reduced. A reduced Nyquist limit increases the risk of aliasing, which results in the misrepre- 
sentation of high velocities. 

[0020] According to a first aspect of the invention, there is provided a method for generating a composite tissue 
image and tissue movement image in an ultrasound system comprising: acquiring echo signals for a plurality of range 
positions along a plurality of ultrasonic beams covering a spatial region of interest during a first time period; acquiring 

50 echo signals for a plurality of range positions along the plurality of ultrasonic beams covering the spatial region of inter- 
est during a second time period; acquiring echo signals for a plurality of range positions along the plurality of ultrasonic 
beams covering the spatial region of interest during a third time period; processing acquired echo signals from at least 
said first and second time periods to produce a first frame of a composite tissue image and tissue movement image; 
and processing acquired echo signals from at least said second and third time periods to produce a second frame of a 

55 composite tissue image and tissue movement image. 

[0021] The method may further comprise: acquiring echo signals for a plurality of range positions along the plurality 
of ultrasonic beams covering the spatial region of interest during a fourth time period; and processing acquired echo 
signals from at least said third and fourth time periods to produce a third frame of a composite tissue image and tissue 
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movement image. 

[0022] Processing acquired echo signals from at least said first and second time periods to produce a first frame of 
a composite tissue image and tissue movement image may include processing acquired echo signals from said third 
time period; and processing acquired echo signals from at least said second and third time periods to produce a second 
5 frame of a composite tissue image and tissue movement strain rate image includes processing acquired echo signals 
from said fourth time period. 

[0023] The tissue movement image may comprise a strain rate image and processing acquired echo signals may 
comprise: estimating strain rates for said range positions inside said spatial region of interest. 
[0024] The step of estimating strain rates may comprise: estimating tissue velocity for range positions along said 
10 ultrasonic beams based on the echo signals; and calculating the strain rate as a spatial derivative of the tissue velocity. 
[0025] The spatial derivative may be found with a linear regression of the tissue velocity for range positions along 
said ultrasonic beams. 

[0026] The step of estimating strain rates may comprise: estimating tissue velocity for range positions along the 
ultrasonic beams based on the echo signals; and calculating the strain rate by determining a velocity difference 
15 between estimated tissue velocities associated with at least a first and second range positions and dividing the velocity 
difference by a distance between the first and second range positions. 

[0027] The step of estimating the strain rate may comprise: estimating a complex pulse-to-pulse correlation R(r) for 
a number of range positions along an ultrasonic beam based on the echo signals; determining a strain correlation func- 
tion, S(r), over a radial distance dr according to an equation S(r)=conj(R(r))*R(r+dr) ; and calculating the strain rate 
20 according to an equation SV(r)= c/(4rcdrTfo) phase (S(r)). 

[0028] The strain correlation function S(r) may be temporally averaged. 

[0029] The step of estimating the strain rate may comprise: estimating a complex pulse-to-pulse correlation for a 
number of range positions along an ultrasonic beam, based on the echo signals; calculating a strain correlation function 
from at least two range positions separated by a given radial distance; and calculating the strain rate based on the 
25 phase of the strain correlation function. 

[0030] The strain correlation function may be given by multiplying the conjugate of the complex pulse-to-pulse cor- 
relation for a first range position by the complex pulse-to-pulse correlation for a second range position where said sec- 
ond range position is located the given radial distance from said first range position. 

[0031] The strain rate may be given by dividing a numerator defined as the product of the phase angle of the strain 
30 correlation function and the speed of sound by a denominator defined as the product of 4, n, the given radial distance, 
the ultrasound frequency and the time between consecutive pulses of said multiple of pulses. 
[0032] The tissue movement image may be a tissue velocity image and the tissue movement image may be a strain 
image. 

[0033] According to a second aspect of the invention, there is provided a method for generating tissue deformation 
35 information comprising: acquiring echo signals for a plurality of range positions along an ultrasonic beam in an area of 
interest to cover a spatial region; estimating a tissue deformation value for said range positions inside said spatial 
region; and displaying tissue deformation values for each range position on a display unit to provide an image of said 
deformation values for said spatial region. 

[0034] The tissue deformation values for each range position may be displayed at spatial coordinates on a display 
40 unit associated with said spatial region, to provide a live, real-time image of said deformation values for said spatial 
region. 

[0035] The tissue deformation value may comprise a strain rate defined as a spatial derivative of tissue velocity. 
[0036] The step of estimating the tissue deformation value may comprise estimating a strain, defined as a strain 
rate integrated over a given time interval. 

45 [0037] The strain rate may be defined as a spatial derivative of tissue velocity. 

[0038] The step of estimating the strain may comprise: estimating a complex pulse-to-pulse correlation for a 
number of range positions along the ultrasonic beam, based on the echo signals; calculating a strain correlation func- 
tion from at least two range positions separated by a given radial distance; calculating the strain rate based on the 
phase of the strain correlation function;, and calculating strain by accumulating the strain rate over the given time inter- 

50 val. 

[0039] The strain correlation function may be given by multiplying the conjugate of the complex pulse-to-pulse cor- 
relation for a first range position by the complex pulse-to-pulse correlation for a second range position where said sec- 
ond range position is located the given radial distance from said first range position. 

[0040] The strain rate may be given by dividing a numerator defined as the product of the phase angle of the strain 
55 correlation function and the speed of sound by a denominator defined as the product of 4, n, the given radial distance 
and the time between consecutive pulses of said multiple of pulses. 
[0041] The given time interval may relate to events in the cardiac cycle. 
[0042] The given time interval may be triggered by an artifact in an ECG trace. 
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[0043] The tissue deformation value may be strain and the step of estimating the tissue deformation value may 
comprise: estimating a strain rate for a given sample volume over a plurality of frames, each of said frames being sep- 
arated by a frame interval; multiplying the estimated strain rate for each of said plurality of frames by the frame interval 
to determine a strain value for each of said plurality of frames; and summing the strain values for each of said plurality 
5 of frames. 

[0044] The tissue deformation value may be strain and the step of estimating the tissue deformation value may 
comprise: estimating a strain rate for a first sample volume over a frame interval; multiplying the strain rate for said first 
sample volume by the frame interval to determine a first strain value; estimating a tissue velocity for said first sample 
volume and calculating a frame-to-frame relative displacement value; estimating a strain rate for a second sample vol- 
w ume over a frame interval, said second sample volume being displaced from said first sample volume by said frame-to- 
frame relative displacement value; multiplying the strain rate for said second sample volume by the frame interval to 
determine a second strain value; and summing at least said first and second strain values. 

[0045] he tissue deformation value may be strain and the step of estimating the tissue deformation value may com- 
prise: estimating a strain rate for a first sample volume over a frame interval; multiplying the strain rate for said first sam- 

15 pie volume by the frame interval to determine a first strain value; estimating a tissue velocity for said first sample volume 
and calculating a frame-to-frame relative displacement value; estimating a strain rate for at least one additional sample 
volume over at least one additional frame interval, said at least one additional sample volume being displaced from said 
first sample volume by said frame-to-frame relative displacement value; multiplying the strain rate for said at least one 
additional sample volume by the frame interval to determine at least one additional strain value; and summing said first 

20 strain value and said at least one additional strain value. 

[0046] According to a third aspect of the invention, there is provided a method for generating tissue deformation 
information for a tissue portion comprising: estimating tissue velocity for a number of sample volumes along an ultra- 
sonic beam, a first end sample volume and a second end sample volume of said sample volumes defining a spatial off- 
set; determining whether said first end sample volume and said second end sample volume are within said tissue 

25 portion; automatically adjusting said spatial offset such that said first end sample volume and said second end sample 
volume are within said tissue portion and said spatial offset may be maximized; and calculating a tissue deformation 
value as a spatial derivative of the tissue velocity over said spatial offset. 

[0047] The step of determining may comprise: determining whether a grayscale value associated with said first end 
sample volume may be above a threshold; and determining whether a grayscale value associated with said second end 
30 sample volume may be above a threshold. 

[0048] The step of determining may comprise: determining whether an absolute power estimate associated with 
said first end sample volume may be above a threshold; and determining whether an absolute power estimate associ- 
ated with said second end sample volume may be above a threshold. 

[0049] The step of determining may comprise: determining whether a magnitude of an auto-correlation function 
35 with unity temporal lag associated with said first end sample volume may be above a threshold; and determining 
whether a magnitude of an auto-correlation function with unity temporal lag associated with said second end sample 
volume may be above a threshold. 

[0050] The step of determining may comprise: determining whether a magnitude of a strain correlation associated 
with said first end sample volume may be above a threshold; and determining whether a magnitude of a strain correla- 

40 tion associated with said second end sample volume may be above a threshold. 

[0051] The tissue deformation value may be strain rate and the strain rate may be used to calculate strain. 
[0052] According to a fourth aspect of the invention, there is provided a method for generating tissue deformation 
information comprising: determining a tissue velocity for a plurality of sample volumes along an ultrasound beam; esti- 
mating a first strain rate as a spatial derivative of tissue velocity based upon a first spatial offset comprising a first subset 

45 of said sample volumes; estimating a second strain rate as a spatial derivative of tissue velocity based upon a second 
spatial, offset comprising a second subset of said sample volumes; and estimating a weighted strain rate based on a 
weighted sum of said first strain rate and said second strain rate. 

[0053] The first strain rate may be weighted in proportion to a strain correlation estimate for said first spatial offset 
and the second strain rate may be proportional to a strain correlation estimate for said second spatial offset. 
so [0054] The first strain correlation estimate may be a function of a first signal correlation estimate and the second 
strain correlation estimate may be a function of a second signal correlation estimate. 

[0055] The method may further comprise: estimating said first and second signal correlation estimates wherein 
said estimation comprises calculating the product of a complex conjugate of a quadrature demodulated Doppler signal 
received from of an end sample volume of said spatial offset during a first temporal sample and a quadrature demodu- 
55 lated Doppler signal received from said end sample volume during a subsequent temporal sample. 

[0056] The method may further comprise: estimating said first and second signal correlation estimates wherein 
said estimation comprises calculating the product of a complex conjugate of a quadrature demodulated Doppler signal 
received from of a first end sample volume of said spatial offset during a first temporal sample and a quadrature demod- 
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ulated Doppler signal received from a second sample volume separated from said first end sample volume by a spatial 
lag during a subsequent temporal sample. 

[0057J The spatial lag may be proportional to the tissue velocity determined for the first end sample volume. 
According to a fifth aspect of the invention, there is provided a method for 

5 [0058] estimating a tissue velocity comprising; transmitting an ultrasound signal at a fundamental frequency; 
receiving an echo of said ultrasound signal from a sample volume and band pass filtering three equal copies of said 
echo, each copy being filtered in a different band pass range producing a first filtered echo, a second filtered echo and 
a third filtered echo respectively, said first filtered echo being band pass filtered to a range centered at a first frequency, 
f h that may be less than a reference frequency related to said fundamental frequency, said second filtered signal being 

10 filtered to a range centered at a second frequency, f 2 , that may be greater than said reference frequency and said third 
filtered signal being filtered to a range centered at a third frequency, f 3 , that may be equal to said reference frequency; 
estimating a difference correlation estimate as the product of a complex conjugate of a signal correlation estimate of 
said first filtered signal and a signal correlation estimate of said second filtered signal; calculating a first tissue velocity 
that may be proportional to the angle of the difference correlation estimate divided by the difference between said sec- 

15 ond and first frequencies; calculating a number of candidate velocities, each of said candidate velocities being propor- 
tional to the sum of the angle of a signal correlation estimate of said third filtered signal and a frequency factor divided 
by the third frequency, where a candidate velocity may be calculated for all values of the frequency factor in a range 
between a negative and positive value of (f 3 - (f 2 -fi))l (2( f 2 - )); and selecting one of said candidate velocities that 
may be closest to said first tissue velocity as an output tissue velocity. 

20 [0059] The reference frequency may be equal to a harmonic frequency of said fundamental frequency. 

[0060] According to a sixth aspect of the invention, there is provided a method for estimating a tissue velocity com- 
prising; transmitting at least two ultrasound pulses along an ultrasound beam; receiving an echo signal for each of said 
transmitted ultrasound signals; band pass filtering each of said received echo signals with at least two different band 
pass filters with center frequencies f 1 and f 2 , producing a first filtered echo signal and a second filtered echo signal; 

25 calculating a pulse-to-pulse correlation estimate for each of said filtered echo signals; and calculating a first tissue 
velocity that may be proportional to the difference in phase angle between the two said correlation estimates, divided 
by the difference between f 2 and f 1 . 

[0061] The method may further comprise: band pass filtering said received echo signals with a third band pass filter 
with center frequency f 3 , producing a third filtered echo signal: calculating a third pulse-to-pulse correlation estimate for 

30 said third filtered echo signal; calculating a number of candidate velocities, each of said candidate velocities being pro- 
portional to the sum of the angle of the correlation estimate for said third filtered echo signal and a frequency factor 
divided by f 3 , where a candidate velocity may be calculated for all values of the frequency factor in a range between a 
negative and positive value of (f 3 - ( f 2 - U ) ' ( 2 ( h ~ f i))> and selecting one of said candidate velocities that may be clos- 
est to said first tissue velocity as an output tissue velocity. 

35 [0062] The third pulse-to-pulse correlation estimate may be calculated with a spatial lag proportional to said first tis- 
sue velocity. 

[0063] At least one of said three band pass filters are centered in a frequency range generated by nonlinear prop- 
agation of said transmitted pulses. 

[0064] According to a seventh aspect of the invention, there is provided a method for performing quantitative stress 
40 echo ultrasound comprising: estimating and storing a tissue deformation value for a heart wall tissue segment of a 

patient over a cardiac interval during each of at least two stress periods, where a level of stress on the patient may be 

different for each of said at least two stress periods; and simultaneously displaying the estimated strain rates for each 

of said at least two stress periods as a function of time over the cardiac interval. 

[0065] The cardiac interval may correspond to an R to R interval of the cardiac cycle. 
45 [0066] The display of strain rates for at least one of said at least two stress periods may be time scaled such that 

the length of the cardiac interval during each of said at least two stress periods appears to be equal in length. 

[0067] The tissue deformation value may be strain rate. 

[0068] The tissue deformation value may be strain accumulated over said cardiac interval. 
[0069] The at least two stress periods may comprise three stress periods. 

50 [0070] According to an eighth aspect of the invention, there is provided a method for generating tissue deformation 
information comprising: acquiring echo signals for a plurality of beams and a plurality of range positions along ultrasonic 
beams in an area of interest to cover a spatial region; determining a beam angle between the ultrasonic beams and a 
principle direction for local tissue deformation; computing at least one angle corrected tissue deformation parameter 
along said principal direction for at least one spatial location; and displaying at least one of the said angle corrected tis- 

55 sue deformation parameters on a display unit. 

[0071] The ultrasonic beams may be generated with a high lateral resolution inside said area of interest. 

[0072] The said beam angle determination may be computed based on a direction along and perpendicular to a 

user defined polygon. 
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[0073] The said computation of at least one angle corrected tissue deformation parameter may comprise: comput- 
ing a radial velocity gradient radially along the ultrasound beam; computing a lateral velocity gradient laterally between 
beams at a fixed range location; and deriving angle corrected tissue deformation parameters as a linear combination of 
said radial and lateral velocity gradients determined by said beam angle. 
5 [0074] The method may further comprise: spatially averaging said radial and lateral velocity gradients. 

[0075] Changes in at least one of said angle corrected tissue deformation parameters may be displayed as a func- 
tion of time for a given anatomical location. 

[0076] The said display may be a M-Mode display displaying at least one of said angle corrected tissue deformation 
parameters with time versus location on said user defined polygon. 

10 [0077] According to a ninth aspect of the invention, there is provided a method for real-time imaging of temporally 
accumulated a tissue movement property comprising: acquiring echo signals for a plurality of ultrasonic beams and a 
plurality of range positions along said ultrasonic beams in an area of interest to cover a spatial region; estimating at 
least one tissue movement property inside the said area of interest; obtaining a sequence of trigger events; accumulat- 
ing said tissue movement property estimates since a most recent trigger event for said spatial region into an accumu- 

15 lation image; and displaying said accumulation image on a display unit in real-time. 

[0078] The tissue movement property may be tissue velocity and said accumulation image may be distance com- 
puted as a time velocity integral since the most recent trigger event. 

[0079] The said tissue movement property may be strain rate and the associated said accumulation image may be 
strain computed as a monotonic map of the sum of strain rate estimates since the most recent trigger event. 

20 [0080] The said trigger events may identify an R event in the cardiac cycle. 

[0081] The method may further comprise: measuring tissue velocity along said ultrasound beams for all range posi- 
tions inside said area of interest; and accumulation tissue velocity from different range positions for each temporally 
successive frame in order to compensate for movement of anatomical locations along the ultrasound beam measured 
by said tissue velocity measurement. 

25 [0082] Thus an ultrasound system and method for calculation and display of tissue deformation parameters are dis- 
closed. 

[0083] In a first aspect of a preferred embodiment of the present invention, an ultrasound acquisition technique that 
allows a high frame rate in tissue velocity imaging or strain rate imaging is disclosed. With this acquisition technique the 
same ultrasound pulses are used for the tissue image and the Doppler based image. A sliding window technique is 
30 used for processing. 

[0084] In a second aspect of an embodiment of the present invention, strain is determined by an accumulation of 
strain rate estimates for consecutive frames over an interval. The interval may be a triggered interval generated by, for 
example, an R-wave in an ECG trace. The strain calculation may be improved by moving the sample volume from which 
the strain rate is accumulated from frame-to-frame according to relative displacement of the tissue within the original 
35 sample volume. The relative displacement of the tissue is determined by the instantaneous tissue velocity of the sample 
volume. 

[0085] In another aspect of a preferred embodiment of the present invention, dr, the spatial offset used in an esti- 
mation of strain rate, is varied adaptively throughout the image. The spatial offset, dr, can be maximized to cover the 
entire tissue segment (e.g., heart wall width) while still keeping both of the sample volumes at each end of the offset 

40 within the tissue segment. This may be accomplished by determining whether various parameters (e.g., grayscale 
value, absolute power estimate, magnitude of the autocorrelation function with unity temporal lag and/or magnitude of 
strain correlation) of the sample volumes within in the spatial offset are above a given threshold. 
[0086] In a still further aspect of a preferred embodiment of the present invention, a generalized strain rate estima- 
tor that is based on a weighted sum of two-sample strain rate estimators with different spatial offsets is employed. The 

45 weights are proportional to the magnitude of the strain rate correlation estimate for each spatial offset, and thus reduce 
the effect of noisy, i.e. poorly correlated, samples. 

[0087] In another aspect of a preferred embodiment of the present invention, an improved signal correlation esti- 
mator that uses a spatial lag in addition to the usual temporal lag is disclosed. The spatial lag is found from the tissue 
velocity. The improved signal correlation estimator can be utilized both in the estimation of strain rate and tissue velocity. 

so [0088] In yet another aspect of a preferred embodiment of the present invention, tissue velocity is estimated in a 
manner that reduces aliasing while maintaining spatial resolution. Three copies of a received ultrasound signal are 
bandpass filtered at three center frequencies. The middle of the three center frequencies is centered at the second har- 
monic of the ultrasound signal. A reference tissue velocity is estimated from the two signals filtered at the outside center 
frequencies. The reference tissue velocity is used to choose a tissue velocity from a number of tissue velocities esti- 

55 mated from the signal centered at the second harmonic. 

[0089] In another aspect of a preferred embodiment of the present invention, a method to estimate the strain rate 
in any direction, not necessarily along the ultrasound beam, based on tissue velocity data from a small region of interest 
around a sample volume is disclosed. 
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[0090] In another aspect of a preferred embodiment of the present invention, a plurality of quantitative tissue defor- 
mation parameters, such as tissue velocity, tissue velocity integrals, strain rate and/or strain, may be presented as func- 
tions of time and/or spatial position for applications such as stress echo. For example, strain rate or strain values for 
three different stress levels may be plotted together with respect to time over a cardiac cycle. Parameters which are 
5 derived from strain rate or strain velocity, such as peak systolic wall thickening percentage, may be plotted with respect 
to various stress levels. 

[0091] The invention will now be described in greater detail, by way of example, with reference to the drawings, in 
which:- 

10 Fig. 1 illustrates the pulse order and beam direction for three different interleave group sizes. 
Fig. 2 illustrates a typical ultrasound acquisition procedure for a tissue Doppler application. 
Fig. 3 illustrates an ultrasound system according to a preferred embodiment of the present invention. 

15 

Fig. 4 illustrates an ultrasound acquisition procedure for a tissue Doppler, strain or strain rate application according 
to a preferred embodiment of the present invention. 

Fig. 5 illustrates a graph of the results of a computer simulated comparison of a linear regression strain rate esti- 
20 mator according to the prior art and a weighted strain rate estimator according to a preferred embodiment of the 
present invention. 

Fig. 6 illustrates a graph of the results of a computer simulated comparison of a linear regression strain rate esti- 
mator according to the prior art and a weighted strain rate estimator according to a preferred embodiment of the 
25 present invention. 

Fig. 7 illustrates the coordinates r, u, v, and w, and the insonation angle a used by a strain rate estimate angle cor- 
rection technique according to a preferred embodiment of the present invention. 

30 Fig. 8 illustrates the velocity components Vy, v^, and v p the distance Ar and the angle a in a small muscle segment 
used by a strain rate estimate angle correction technique according to a preferred embodiment of the present 
invention. 

Fig. 9 illustrates a display of strain rate from multiple stress levels as a function of time for a normal case according 
35 to a preferred embodiment of the present invention. 

Fig. 1 0 illustrates a display of accumulated strain from multiple stress levels as a function of time for a normal case 
according to a preferred embodiment of the present invention. 

40 Fig. 1 1 illustrates a display of strain rate from multiple stress levels as a function of time for a pathological case 
according to a preferred embodiment of the present invention. 

Fig. 12 illustrates a display of accumulated strain from multiple stress levels as a function of time for a pathological 
case according to a preferred embodiment of the present invention. 

45 

Fig. 13 illustrates a display of a strain derived parameter, peak systolic wall thickening, as a function of stress level 
according to a preferred embodiment of the present invention. 

[0092] A method and apparatus are described for generating diagnostic images of tissue deformation parameters, 
so such as strain rate, strain and tissue velocity, in real time and/or in a post-processing mode. In the following description, 
numerous specific details are set forth in order to provide a thorough understanding of the preferred embodiments of 
the present invention. It will be apparent, however, to one of ordinary skill in the art that the present invention may be 
practiced without these specific details. 

[0093] A block diagram for an ultrasound imaging system according to a preferred embodiment of the present 
55 invention is shown in Fig. 3. A transmitter 140 drives an ultrasonic transducer 142 to emit a pulsed ultrasonic beam 144 
into the body. The ultrasonic pulses are back-scattered from structures in the body, like muscular tissue, to produce ech- 
oes which return to and are detected by the transducer 142. A receiver 146 detects the echoes. The echoes are passed 
from the receiver 146 to a complex demodulation stage 148 and a tissue processing stage 149. The complex demodu- 
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lation stage 148 demodulates the echo signals to form I, Q data pairs representative of echo signals. The demodulated 
I, Q data pairs are complex Doppler signals that are passed to a tissue deformation calculation stage 150 which carries 
out tissue velocity, strain rate and/or strain calculations as explained below. The complex Doppler signal is associated 
with a sample volume defined by a range position and beam in a region of interest. A complex Doppler signal typically 
5 may comprise a segment of data samples which is used to estimate the Doppler shift. The echo signals are also passed 
to the tissue processing stage 1 49, which performs processing such as B-mode processing to form a 2D or 3D image 
of the anatomical structure scanned. 

[0094] The tissue deformation values, e.g., tissue velocity, strain rate and/or strain, output by the tissue deformation 
calculation stage 150 and the tissue image values output by the tissue processing stage 149 are passed to a display 
w system 152 for display. The display system 152 includes a monitor 154. 

[0095] United States Patent Application No. 09/167,896, filed October 7, 1998 and entitled "A METHOD AND 
APPARATUS FOR PROVIDING REAL-TIME CALCULATION AND DISPLAY OF STRAIN IN ULTRASOUND IMAG- 
ING," describes a manner in which a strain rate may be estimated using the system of Fig. 3. 

[0096] For Strain Rate Imaging (SRI) and other Doppler based applications where a low pulse repetition frequency 
15 (PRF) is acceptable, a scanning procedure that allows higher frame rate may be used. Instead of collecting separate 
tissue frames as illustrated in Fig. 2, the number of beams in the Doppler subframes can be increased to allow tissue 
visualization based on only these frames. The acquisition of separate tissue frames becomes unnecessary. Fig. 4 illus- 
trates a scanning procedure that allows a high frame rate. This scanning procedure may be used in either tissue Dop- 
pler or SRI applications. In the example of Fig. 4, the packet size N = 3 and the interleave group size N jnt = N b . T is 
20 the pulse repetition time, t T and t D are the times needed to acquire a tissue frame and a Doppler frame respectively, and 
t F is the total acquisition time for one tissue Doppler or SRI frame. As illustrated in Fig. 4, a Doppler frame is still gen- 
erated from N subframes (the subframes are numbered 160, 161, 162, 163 and 164), but a sliding window technique 
may be used, so the time to produce one Doppler or SRI frame will be only 



25 
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55 



t FSRJ ~~ *T> 



((8) 



30 assuming that the time to acquire one Doppler subframe is equal to the time to acquire one tissue frame in the conven- 
tional method. Comparing equations (( 7 ) and (( 8 ) one can see that the acquisition time for one frame is greatly 
reduced and, thus, allowing a higher frame rate. 

[0097] One parameter that may be calculated by the tissue deformation calculation stage 1 50 is strain. The relation 
between the strain and the strain rate can be developed by way of an example. Consider a one-dimensional homoge- 
35 neous object of length L(t) that has a spatially constant strain rate field s(t). The term "strain rate" is here used for the 
velocity gradient. The velocity field is thus given as: 



v(/,r) = s(r)r, 



((9) 



where r is the position in the object. The velocity at r = 0 is set to zero for simplicity, but the same relations will apply 
45 also when v(t,0) differs from zero. 

[0098] The change in length over a small time step Af can then be estimated as 



I(/ + A/)-Z(0«A/j(/)iL(0. 



[0099] Letting At ^ Owe get the temporal derivative of the length: 



((10) 
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dt ^ A/ 



((11) 



70 [0100] The solution to this differential equation is 



I(/) = I 0 expfcj(r)<*r) 



75 



((12) 



and the strain is finally found as 



20 



25 



e (t) = L ° . 100% = [exp(£ 5(r)rfr)- 1]« 100%. 



((13) 



[0101] The strain e(i) in a sample volume in the image can be estimated in real-time by replacing the integration in 
equation (( 13 ) with a cumulative sum: 



30 



35 



e{i) = [exp(C(/))-l]l00%, 
C(0 = C(i-l) + j(i)A/. 



((14) 



40 



45 



[0102] Here / is the frame number and At is the time between each frame. C(i) is the cumulative sum, and s(i) is 
the strain rate estimate for the given sample volume. The accumulation can also be reset at any time, for instance at a 
specific time trigged by an ECG-signal, by setting C(i-1) to zero for the corresponding frame number /. The calculation 
above can be performed for every sample volume in the image, and the visualization can be performed in the same way 
as for tissue velocity imaging (TVI) and SRI, only using a color map representing strain rather than tissue velocity or 
strain rate. 

[0103] A further improvement is possible if the tissue velocity v is also available for each sample volume. In the 
cumulative sum for radial sample volume number m 0 , the strain rate estimate might then be taken from a different sam- 
ple volume given by the tissue velocity. First, the frame-to-frame relative dis 



50 



d = vAtk t 



((15) 



55 placement index is estimated as where v is the tissue velocity estimate in sample number m 0 , and k s is the spatial sam- 
pling frequency. Next, the strain rate estimate from the sample volume number 
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m = m Q +d 



((16) 



5 



is used in the cumulative sum, rather than m 0 . If the tissue movement is only in the direction of the beam, this method 
allows the cumulative summation to track the motion of the same anatomical sample during its movement. Even if the 
tissue movement is in other directions, an improvement is expected. 
10 [0104] The strain rate estimator in Application No. 09/167,896 was in its simplest form described as: 



where r is the radial position along an ultrasound beam, v is the tissue velocity, and dr is the spatial offset. This spatial 
offset can be varied adaptively throughout the image. Given an upper and lower limit on the size of dr, it can be 
increased as much as possible while still keeping both of the sample volumes at each end of the offset within the tissue. 

20 There are several different criteria that can be used to ensure that the offset is within the tissue. One possible criteria is 
that the corresponding tissue sample volumes must have a grayscale value above a given limit. Another possible crite- 
ria is that the power estimates of the sample volumes must have absolute values above a given limit. Another possible 
criteria is that in either of the two sample volumes, the magnitude of the auto-correlation function with unity temporal lag 
must be above a given limit. Another possible criteria is that the magnitude of the strain correlation (described in equa- 

25 tion (8) in Application No. 09/167,896) must be above a given limit. Any of these criteria one can be used separately, or 
they can be combined so that two or more criteria must be met for a positive determination that the sample volumes at 
the end of the offset dr are within the tissue. 

[0105] The tissue deformation calculation stage 14 may calculate strain rate using a strain rate estimator that is 
based on several samples, and is weighted with the magnitude of a strain correlation estimate. Consider a quadrature 
30 demodulated Doppler signal x(m,n), where m is the spatial sample volume index, and n is the temporal index. The sig- 
nal is assumed to have been acquired using a center frequency f 0 , a pulse repetition time 7, and a radial sampling fre- 
quency r s equal to the radial size of the point spread function. The speed of sound in the imaged object is assumed to 
be c. An estimator for strain rate based on M spatial and N temporal samples of the Doppler signal is: 



j(r) = (v(r + </r)-v(>) )ldr 



((17) 



15 



35 




((18) 



40 



45 where 




((19) 



50 



is the strain rate correlation estimate, 



55 
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<0,(m)=-ZS(m), 



((20) 



5 



m 



is the angle of the strain rate correlation estimate, and 



10 



2 




((21) 



m 



15 



is a weighing factor. The signal correlation estimate R(m) is described below. 

[0106] The strain rate estimator of equation (18) has certain advantages over the prior art Myocardial Velocity Gra- 
20 dient (MVG) technique first described in D. Fleming et al., "Myocardial velocity gradients detected by Doppler imaging" 
Br. J. Radiol., 67(799):679-688, 1994, and further developed by Uematsu et al., "Myocardial velocity gradient as a new 
indicator of regional left ventricular contraction: Detection by a two-dimensional tissue Doppler imaging technique" J. 
Am. Col. Cardiol., 26(1):217-23, 1995. For example, Fleming and Uematsu disclose the use of a least squares linear 
regression of the velocity data to obtain the velocity gradient (strain rate). Linear regression puts equal weight to all the 
25 velocity samples. The weighted strain estimator of equation (18), however, uses weights that vary with the magnitude 
of the strain rate correlation of equation (19) resulting in an improved strain rate estimation. 

[0107] Figs. 5 and 6 illustrate a computer simulation comparison of a least squares linear regression estimator and 
the strain rate estimator of equation (1 8). Fig. 5 illustrates a linear regression fit (dashed line) and a weighted strain rate 
linear fit (solid line) for simulated velocity estimates (circles) at varying depths. Signals including noise were generated 

30 with a velocity gradient (strain rate) of 1.0 s" 1 . A typical outcome is presented in Fig. 5. Note that the two outermost 
points give a large error for the linear regression line (dashed line), while the effect on the weighted strain rate estimator 
is much less. In Fig. 6, strain rates estimated with the linear regression method (stars) and with the weighted strain 
weight estimator (circles) are compared for 50 independent simulations. Once again, signals including noise were gen- 
erated with a velocity gradient (strain rate) of 1.0 s" 1 . The weighted strain rate estimator shows less variance than the 

35 linear regression method. 

[0108] The signal correlation R (m) (used in equation (1 9) above) can be estimated in different ways. For example, 
one estimate is 



[0109] Spatial averaging may also be used to reduce the variance of R (m) in equation (22) and other estimators 
of R (m)described herein. 

[0110] A more robust method to estimate the signal correlation R(m) is to introduce a spatial lag Am, and correlate 
so signal samples from not just the same depth m, but also from m + Am: 



40 




((22) 



n-1 



45 




((23) 
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[0111] The spatial lag Am preferably is chosen to maximize the magnitude of R(m). One way to chose a Am is 
through a phase root seeking technique such as described in A. Peasvento and H. Ermert, "Time-efficient and exact 
algorithms for adaptive temporal stretching and 2D-correlation for elastographic imaging using phase information" Proc. 
5 of the 1998 Ultrasonic Symposium, to be published, 1998. Alternatively, the inventors have found that the peak magni- 
tude of R(m) is found when the spatial lag Am is chosen equal to the translation of the tissue from pulse to pulse: 



10 



Am = 



PRF 



((24) 



15 where v is the tissue velocity, PRF is the pulse repetition frequency and k s is the spatial sampling frequency of the sig- 
nal. This method requires that an unaliased velocity estimate is available. 

[0112] The tissue deformation calculation stage 1 50 may calculate a velocity estimate as follows. Three equal cop- 
ies of the received signal are band pass filtered with three different fitters. Two narrow band filters centered at f 1 and 
and a third wider band filter centered at f 3 are used, where f 1 <f 3 < f 2 , and f 3 is centered around the second harmonic 
20 component of the signal. The signal correlation of each of these three signals are estimated using equation ( 22 ), 
resulting in the correlation estimates &i(m), R2(m) and ^(m), respectively. 
[01 13] The tissue velocity can be found from the angle of R 3 (m) as: 



25 



30 



. cPRF yt t , 



((25) 



where c is the speed of sound. Unfortunately, the velocity estimate of equation (25) is easily aliased. A difference cor- 
relation is next found as 



35 



k 4 {m) = J*,'(m)* 2 (m). 



((26) 



40 

[0114] The velocity of the tissue is found from the angle of this difference correlation as 



cPRF , 6 . 

as = Tin- — 7\ ("*)» 



((27) 



50 where c is the speed of sound. This velocity estimate is not as easily aliased since (f 2 -fj)< h- However, with equation 
(27) the spatial resolution is poor since narrow band signals were used in the estimation of R -j(m) and R2( m )- To this 
point, this two-frequency velocity estimation method is similar to a method described in Dousse etal., "Two years expe- 
rience in measuring velocities beyond the Nyquist limit with Color Flow Mapper" Proceedings of EURODOP'92, page 
219, Brighton, United Kingdom, 1992. 

55 [0115] To regain the spatial resolution of the estimate in equation (25), the following algorithm is used: For each 
(possibly aliased) velocity estimate £ 3 , several candidate velocities are found as 
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cPRF 



(z£ 3 (m) + 2Jbr)l - 




((28) 




[01 16] Next, the candidate velocity v 3 k that is closest to the (unaliased) difference velocity estimate v d is chosen 
as the output velocity estimate. This way, the spatial resolution of the v 3 estimates is kept, while avoiding the problem 
of aliasing. 

[01 17] A method for angle correction of a strain rate estimation is described with respect to Fig. 7. Locally for each 
muscle segment of the left cardiac ventricle, the coordinates are defined as: 

r - along the ultrasound beam, positive away from the transducer 

I - lateral (beam-to-beam), positive from left to right in the ultrasound image 

u - circumferential, clockwise seen from the apex 

v - meridional (longitudinal), from apex to base 

w - transmural, from endo- to epi-card 

where u, v and w will be approximately perpendicular, as shown in Fig. 7. The strain rates in these directions are termed 
*r> s b s u> s * and s w» respectively. The origin (u,v,w)=(0,0,0) does not need to be defined in relation to the macroscopic 
ventricle geometry, and can be chosen anywhere in the imaged muscle segment. 

[01 18] Furthermore, a is defined as the angle between the v-axis and the r-axis, so that zero degrees corresponds 
to measuring along the muscle in the meridional direction. It is assumed that the angle a is in the v-w-plane (long axis 
or apical views), so the problem becomes two-dimensional. Note that the angle a is negative in Fig. 7. 
[0119] Without loosing generality, it can be assumed that the point (v,w)=(0,0) is not moving. If the strain rate is spa- 
tially homogeneous over a small distance Ar in the muscle segment, the muscle point (v,w) will then move with the 
velocity components: 



[0120] These velocity components are shown in Fig. 8. Fig. 8 is an illustration of the velocity components Vy, v w and 
v r , the distance Ar and the angle a in a small muscle segment. All the parameters are drawn positive, but notice that the 
angle a is usually negative when imaging from the apex, and that Vy, and consequently v r , normally are negative during 
systole. The lateral (beam-to-beam) l-axis is also included for reference. The velocity component along the ultrasound 
beam in position (v,w) becomes 



((29) 



and 



((30) 



K = 



vs v cosa + ws w sin a. 



((31) 
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[0121] Notice that the velocity v r for simplicity is defined positive away from the transducer, i.e., in positive redirec- 
tion. This is opposite of the usual definition for the velocity sign in Doppler imaging. 

[0122] By using velocity-information from more than one beam at a time, it is possible to calculate the strain rate in 
other directions than along the beam. The beams are assumed to be parallel in the region of interest. The vw-axis sys- 
tem is then a rotation of the Ir-axis system by an angle of (a-7i/2), and one can write 



10 



v = rcosa + lsinac 
w = rsina-lcosa 



((32) 



15 



[0123] Substituting these equations in equation (31 ) one gets 



20 



v r = s y (r cos a + 1 sin a) cos a + s w (r sin a - 1 cos a) sin a 



((33) 



25 [0124] Taking the derivatives in the two directions r and I, one gets the two equations 



30 



35 



dr 

dv, 



= s„ cos 2 a + s„ sin 1 a 



r _ 



a 



= j„ sin a cos a -$„ sin a cos a 



((34) 



[01 25] Solving for s v and s w gives 
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*'-* + -a , - B 

dv, dv. 

s w = — '- L cota 

w dr dl 



((35) 
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[0126] This means the strain rates in the anatomical directions v (meridional) and w (transmural) can be found from 
the radial and lateral gradients of the measured radial velocity, as long as the angle a is known. The image plane Ir must 
coincide with the vw plane, which is the case for all apical views and the parasternal long axis view (PLAX). Notice that 
when imaging from the apex, the angle a will be close to zero for most of the ventricle. 

[0127] The same formulas apply if one substitutes v with u, so the strain rate in the u-direction (circumferential) can 
also be found. The image plane Ir must then coincide with the uw plane, which is the case for the short axis view (SAX). 
[0128] There will be some angles where the strain rates are unavailable, though. For the u or the v directions these 
are the angles where tan a approaches infinite values. For the w-direction these are the angles where cot a approaches 
infinite values. 

[0129] In the SAX view and using a sector scan, an approximation of a can automatically be found if the user 
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defines the center of the ventricle. By assuming that the SAX cross section of the ventricle is circular, a at a particular 
location is given as 



a = — -9 b +0 c 



((36) 



10 

where 8^ is the angle of the ultrasound beam that intersects the point (6 b = 0 is defined as the center beam), and 6 C 
is the angle between the center ultrasound beam an imaginary beam from the center of the ventricle through the point. 
[0130] A preliminary test has been performed using this two-dimensional angle correction method. A velocity data 
set from a healthy volunteer was obtained using tissue Doppler imaging with high beam density. The short axis view 

15 was used and the circumferential and transmural strain rate components in three phases of the cardiac cycle (mid sys- 
tole, early diastolic relaxation and mid diastole) were estimated. The myocardium was segmented manually. As 
expected, the resulting images showed that the radial strain rate is equal to the transmural strain rate at twelve and six 
o'clock, and the circumferential strain rate at two and ten o'clock. Except from where the cot a or tan a approach infinity, 
the apparent noise in the images did not seem to be increased by this procedure. 

20 [0131] It is also possible to perform three-dimensional angle correction. Locally for each muscle segment of the left 
cardiac ventricle, the coordinates are defined as: 

x - azimuthal (perpendicular to the image plane) 

25 y - lateral (beam-to-beam) 

z - along the ultrasound beam 

u - circumferential, clockwise seen from the apex 

30 

v - meridional (longitudinal), from apex to base 



w - transmural, from endo- to epi-card 



35 where the triplets x, y, z and u, v, w locally are assumed to be perpendicular. The strain rates in these directions are 
termed s u , s v and respectively. The origin (u,v,w)=(0,0,0) does not need to be defined in relation to the macroscopic 
ventricle geometry, and can be chosen anywhere in the imaged muscle segment. 

[0132] Without loosing generality, it is assume that the point (u,v,w)=(0,0,0) is not moving. If the strain rate is spa- 
tially homogeneous over a small distance Ar in the muscle segment, the muscle point (u.v.w) will then move with the 
40 velocity components: 



v u =us u , v v = vj v , V W =W5 W . 



((37) 



45 



[0133] Using velocity-information from more than one beam at the time it is possible to calculate the strain rate in 
other directions than along the beam. The beams are assumed to be parallel in the region of interest. 
so [0134] Based on formulas for axis rotation it is possible to express the velocity components in the xyz-directions 
rather than in the uvw-directions. The velocity component in the z-direction (along the ultrasound beam), v z , is the one 
found using Tissue Velocity Imaging. The gradients of this velocity component in each of the three spatial directions are 
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v zr=~, r = x,y,z 
or 



((38) 



[0135] The relation to the strain rates in the uvw-directions is 



10 



15 



= A(a,fi,r 



((39) 



20 where A(a,p,v) is a matrix that describes the 3D axis rotation between the uvw-system and the xyz-system, and a, p\ 
and y are the rotation angles about the u-, v- and w-axis respectively. Except for certain rotation angles, this matrix can 
be inverted, and the strain rates can be found as: 
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((40) 



= A- , {a,fi, r ) v, 



30 



[0136] Estimates for the strain rates in the uvw-directions are found by inserting velocity gradient estimators based 
on recorded tissue velocity data. Examples of estimators for the velocity gradients are 



35 



_ v,(r + Ar)-v,(r) 
Ar 



40 



((41) 



where Ax, Ay, and Az are the sampling distances in the azimuth, lateral and radial directions in the ultrasound data 
respectively. Similar methods as described for 1D strain rate can also be used to estimate these velocity gradients, 
45 where the radial increment is replaced by an increment in the x-and y-directions. 

[0137] A further improvement can be achieved by performing a least squares fit of the estimated strain rates to the 
incompressibility equation 



50 



S„ + S„ + S„ ss 0 



((42) 



55 since muscular tissue can be considered incompressible. 

[0138] In two dimensions the strain rate estimates reduce to 
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*u =V a +V^COty0 

= v„+v v tan£ 



((43) 



5 



10 



for images in the uw-plane (short axis images) and 



s v =v a +v ty cota 
s w =v a + Viy tana 



((44) 
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20 for images in the vw-plane (apical images). There will be some angles where the strain rates are unavailable, though. 
For the u or the v directions these are the angles where tan approaches infinite values. For the w-direction these are 
the angles where cot approaches infinite values. 

[0139] The tissue deformation calculations described herein are suited for quantitative stress echo applications. 
There are at least four main quantitative parameters that may be extracted, including: tissue velocity, which quantifies 

25 wall motion; tissue velocity time integrals which quantify accumulated wall motion during a time interval such as systole; 
strain rate (velocity gradient), which quantifies the local wall thickening at a given time instant; and strain (integrated 
strain rate) which quantifies local wall thickening during a time interval such as systole. These parameters are functions 
of both spatial position and time. From these parameters, other clinically relevant parameters may be derived. One way 
to present these parameters is to plot pairs of the parameters against each other (similar to pressure-volume-loops). 

30 Another useful way to present these parameters is to estimate and record (in a cine loop for example) one or more 
parameters from different stress levels in the stress test and then display the respective parameters from the varying 
stress levels simultaneously. 

[0140] During a stress echo examination one of the crucial things to assess is segmental wall motion. Typically, the 
left ventricle is subdivided into segmental areas and a visual assessment of wall motion is done in each of these seg- 

35 ments from the various cine loops that are acquired. The 16 segment ASE model of the left ventricle is currently the 
most common way to subdivide the left ventricle for scoring of stress echo exams. In a visual assessment a given seg- 
ment is compared in terms of motion and wall thickening by visual comparison of similar views (2-chamber, 4-chamber, 
lax, sax, etc.) at different stress levels. The stress levels usually include rest, 1 or more levels of stress induced by exer- 
cise or pharmacological infusion and finally recovery. A normal reading of a segment is that both wall motion and local 

40 wall thickening increase during systole as a function of the applied stress level. 

[0141] Fig. 9 illustrates how time traces of strain rate for a given location or wall segment can be combined from 
multiple stress levels. Strain rates estimated during rest (line 200), medium stress (line 202) and peak stress (line 204) 
are plotted with respect to time. An ECG trace (line 206) is provided for reference at the bottom of the display. A differ- 
ence in heart rate from the various stress levels is accounted for in this example by time scaling the different strain rate 

45 traces. This combined display contains useful clinical information about the local wall function and how the wall seg- 
ment responds to an increase in stress level. The example is a typical normal reading of longitudinal shortening that 
can be recorded with an apical view. It should be noted that the longitudinal shortening assessed in this manner also 
indirectly describes wall thickening in short axis views because of conservation of mass and incompressibility of myo- 
card. The example illustrates a normal reading where longitudinal shortening increases with the stress level. 

so [0142] Fig. 10 is identical to Fig. 9 except that accumulated strain is plotted for rest (line 210), medium stress (line 
212) and peak stress (line 214) instead of the instantaneous strain rate. The Fig. 10 demonstrates how the longitudinal 
shortening increases as a function of stress level. Figs. 11 and 12 correspond to Figs. 9 and 10, respectively, except 
that Fig. 1 1 illustrates a typical pathological reading of strain rates for rest (line 230), medium stress (line 232) and peak 
stress (line 234), and Fig. 12 illustrates a typical pathological reading of accumulated strain for rest (line 240), medium 

55 stress (line 242) and peak stress (line 244). The example of Figs. 11 and 12 illustrates a case with normal rest values 
for longitudinal shortening, but with a reduction in shortening when the stress level increases. At peak stress the curves 
illustrate a reverse in both strain rate and strain which can indicate passive stretching of the local wall segment. 
[0143] Fig. 13 illustrates how characteristic values extracted from the strain information for a given location or wall 
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segment can be displayed as a function of stress level. The example in Fig. 1 3 is the maximal systolic longitudinal short- 
ening which is plotted as a function of stress level. The normal case (line 250) with a uniform increase in longitudinal 
shortening and the pathological case (line 252) with a decrease in longitudinal shortening and even a switch to passive 
stretching during systole are illustrated. 
5 [0144] Another useful way to present the quantitative parameters derived from strain is in a Bulls-eye plot by either 
numerically or graphically labeling each of the areas corresponding to LV segments according to the associated strain 
derived parameters. The values illustrated in Fig. 13 are examples of such useful strain derived parameters. 

Claims 

10 

1. A method for generating a composite tissue image and tissue movement image in an ultrasound system compris- 
ing: 

acquiring (146) echo signals for a plurality of range positions along a plurality of ultrasonic beams (144) cover- 
15 ing a spatial region of interest during a first time period; 

acquiring (146) echo signals for a plurality of range positions along the plurality of ultrasonic beams (144) cov- 
ering the spatial region of interest during a second time period; 

20 acquiring (146) echo signals for a plurality of range positions along the plurality of ultrasonic beams (144) cov- 

ering the spatial region of interest during a third time period; 

processing (149, 150) acquired echo signals from at least said first and second time periods to produce a first 
frame of a composite tissue image and tissue movement image; and 

25 

processing (149, 1 50) acquired echo signals from at least said second and third time periods to produce a sec- 
ond frame of a composite tissue image and tissue movement image. 

2. The method according to claim 1 , further comprising: 

30 

acquiring (146) echo signals for a plurality of range positions along the plurality of ultrasonic beams (144) cov- 
ering the spatial region of interest during a fourth time period; and 

processing (149, 150) acquired echo signals from at least said third and fourth time periods to produce a third 
35 frame of a composite tissue image and tissue movement image. 

3. A method for generating tissue deformation information comprising: 

acquiring (146) echo signals for a plurality of range positions along an ultrasonic beam (144) in an area of inter- 
40 est to cover a spatial region; 

estimating (150) a tissue deformation value for said range positions inside said spatial region; and 

displaying (152) tissue deformation values for each range position on a display unit to provide an image of said 
45 deformation values for said spatial region. 

4. A method for generating tissue deformation information for a tissue portion comprising: 

estimating (150) tissue velocity for a number of sample volumes along an ultrasonic beam (144), a first end 
so sample volume and a second end sample volume of said sample volumes defining a spatial offset; 

determining whether said first end sample volume and said second end sample volume are within said tissue 
portion; 

55 automatically adjusting said spatial offset such that said first end sample volume and said second end sample 

volume are within said tissue portion and said spatial offset is maximized; and 

calculating a tissue deformation value as a spatial derivative of the tissue velocity over said spatial offset. 
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A method for generating tissue deformation information comprising: 

determining a tissue velocity for a plurality of sample volumes along an ultrasound beam (144); 

estimating (150) a first strain rate as a spatial derivative of tissue velocity based upon a first spatial offset com- 
prising a first subset of said sample volumes; 

estimating a second strain rate as a spatial derivative of tissue velocity based upon a second spatial offset 
comprising a second subset of said sample volumes; 

estimating a weighted strain rate based on a weighted sum of said first strain rate and said second strain rate. 
A method for estimating a tissue velocity comprising; 

transmitting (140) an ultrasound signal at a fundamental frequency; 

receiving (146) an echo of said ultrasound signal from a sample volume and band pass filtering three equal 
copies of said echo, each copy being filtered in a different band pass range producing a first filtered echo, a 
second filtered echo and a third filtered echo respectively, said first filtered echo being band pass filtered to a 
range centered at a first frequency, f 1f that is less than a reference frequency related to said fundamental fre- 
quency, said second filtered signal being filtered to a range centered at a second frequency, f 2 , that is greater 
than said reference frequency and said third filtered signal being filtered to a range centered at a third fre- 
quency, f 3 , that is equal to said reference frequency; 

estimating a difference correlation estimate as the product of a complex conjugate of a signal correlation esti- 
mate of said first filtered signal and a signal correlation estimate of said second filtered signal; 

calculating a first tissue velocity that is proportional to the angle of the difference correlation estimate divided 
by the difference between said second and first frequencies; 

calculating a number of candidate velocities, each of said candidate velocities being proportional to the sum of 
the angle of a signal correlation estimate of said third filtered signal and a frequency factor divided by the third 
frequency, where a candidate velocity is calculated for all values of the frequency factor in a range between a 
negative and positive value of (f 3 - (f 2 - /?)) / (2( f 2 - )); and 

selecting one of said candidate velocities that is closest to said first tissue velocity as an output tissue velocity. 
A method for estimating a tissue velocity comprising; 

transmitting (140) at least two ultrasound pulses along an ultrasound beam (144); 
receiving (146) an echo signal for each of said transmitted ultrasound signals; 

bandpass filtering each of said received echo signals with at least two different bandpass filters with center fre- 
quencies f 1 and f 2 , producing a first filtered echo signal and a second filtered echo signal; 

calculating a pulse-to-pulse correlation estimate for each of said filtered echo signals; and 

calculating a first tissue velocity that is proportional to the difference in phase angle between the two said cor- 
relation estimates, divided by the difference between f 2 and ff. 

A method for performing quantitative stress echo ultrasound comprising: 

estimating (150) and storing a tissue deformation value for a heart wall tissue segment of a patient over a car- 
diac interval during each of at least two stress periods, where a level of stress on the patient is different for each 
of said at least two stress periods; and 

simultaneously displaying (152) the estimated strain rates for each of said at least two stress periods (200, 
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202) as a function of time over the cardiac interval. 

9. A method for generating tissue deformation information comprising: 

5 acquiring (146) echo signals for a plurality of beams and a plurality of range positions along ultrasonic beams 

(144) in an area of interest to cover a spatial region; 

determining a beam angle between the ultrasonic beams and a principle direction for local tissue deformation; 

w computing (1 50) at least one angle corrected tissue deformation parameter along said principal direction for at 

least one spatial location; and 

displaying (152) at least one of the said angle corrected tissue deformation parameters on a display unit (154). 

15 10. A method for real-time imaging of temporally accumulated a tissue movement property comprising: 

acquiring (146) echo signals for a plurality of ultrasonic beams (144) and a plurality of range positions along 
said ultrasonic beams in an area of interest to cover a spatial region; 

20 estimating (1 50) at least one tissue movement property inside the said area of interest; 

obtaining a sequence of trigger events; 

accumulating said tissue movement property estimates since a most recent trigger event for said spatial region 
25 into an accumulation image; and 

displaying (152) said accumulation image on a display unit (154) in real-time. 
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FIG. 13 
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